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Technical Report for NASA grant C036846 (NAG3-2134), 

R.G. Larson, P.I., July 25, 2002 

In this report, we summarize work carried out under the above-named grant, 
primarily by post-doc Hua Hu, and partially by grad students Lei Li and Manish Chopra. 
The work includes studies on droplet evaporation and its effects on temperature and 
velocity fields in an evaporating droplet, new 3-D microscopic particle image 
velocimetry and direct visualization on wall slip in a surfactant solution. With the 
exception of the slip measurements, these projects were those proposed in the grant 
application. Instead of slip flow, the original grant proposed imaging electro-osmotic 
flows. However, shortly after the grant was issued, the PI became aware of work on 
electro-osmotic flows by the group of Saville in Princeton that was similar to that 
proposed, and we therefore elected to carry out work on imaging slip flows rather than 
electro-osmotic flows. The following paragraphs are the detailed descriptions on these 
projects. 

1. Droplet evaporation 

The evaporation of a sessile droplet with a pinned contact line was investigated by 
experiments, by analytic theory, and by computation using the finite element method 
(FEM). We compared the results obtained by our FEM analysis with an analytical 
solution and derived a very simple approximate evaporation rate expression 
m(t) = -n RD( 1 - H)c v ( 0.270 1 2 + 1 .30) for any initial contact angle _ between 0 and _/2 

with in radians. The approximate expression was also compared with droplet 
evaporation data from the literature, and good agreement found without any parameter 
fitting. It was found both theoretically and experimentally that the net evaporation rate 
from the droplet remained almost constant with time for a small initial contact angle (0< 
40°), even though the evaporation flux becomes more strongly singular at the edge of the 
droplet as the contact angle decreases during evaporation. The results of this study were 
applied almost immediately by the group of Kate Stebe at Johns Hopkins University, who 
found that the above formula applies even when evaporation is strongly suppressed by a 


surfactant-monoloayer covering the droplet, as long as the droplet does not become too 
small 


Key Result : A simple and accurate formula for predicting the contact-angle 
dependence of the evaporation rate from a sessile droplet, which has been used by Prof. 
Stebe at Johns Hopkins in research on monolayer-covered droplets. 

2. Temperature and velocity fields in the evaporating droplet 

With the results of the droplet evaporation, we used FEM analysis to compute the 
temperature field in the evaporating droplet. We found that the temperature field varies 
with the contact angle and the surface temperature becomes non-uniform along the 
droplet surface. We believe a surface-tension gradient generated by this non-uniform 
temperature distribution produces a Marangoni flow in the droplet. 

Thus, the 3-D time-dependent flow field produced by the evaporating droplet was 
studied theoretically and experimentally. We developed an analytic lubrication theory for 
this flow in the evaporating droplet both with and without Marangoni stresses. Since the 
most available commercial software doesn’t allow us to solve a moving boundary 
problem with Marangoni stress, we developed our own FEM code to solve the Stokes 
equation in the evaporating droplet. The FEM results confirmed the validity and accuracy 
of our approximate analytic solutions for flow field in the droplet. We found that large 
thermally-induced Marangoni flows are predicted for evaporation of pure water, but that 
these Marangoni flows are severely attenuated in the experiments, evidently due to trace 
contaminants producing an offsetting Marangoni stress due to surface-active agents. 
Defocused particle-tracking velocimetry is presented and applied to map the 3- 
dimensional time-dependent particle displacements. The measured velocity field shows 
that there is a weak recirculation near the droplet surface, evidently due to Marangoni 
stresses produced by temperature gradients arising from evaporation. 

Key Result: An analytic lubrication analysis of the entire 3D flow field in an 
evaporating droplet, including the effects of Marangoni flow, with verification of the 
accuracy of the analytic analysis using full FEM, including effects of evaporative 



cooling. This result was recently used to predict DNA deposition from an evaporating 
droplet, a technique used by the David Schwartz group in Wisconsin for genomic 
mapping, including mapping the genome of the parasite plasmodium that causes malaria. 

3. Microscopic 3-D PIV 

In order to map the micro-fluidic flow, we developed two microscopic versions of 
particle image velocimetry. For the flow in the evaporating droplet, 3-D microscopic 
particle-tracking technique was established. In our experiments, we observed particle ring 
patterns, from which we can extract the vertical particle displacement. By tracking the 
particle, we were able for the first time to obtain the 3-dimesional-velocity field from a 2- 
dimensional particle image. Considering the differences between microscopic and 
macroscopic flow properties, we built our own microscopic PIV(particle-image- 
velocimetry) for general 2-D micro-fluidic system, which we then employed to study the 
wall-slip velocity in a surfactant solution. 

Key Result: A new particle imaging velocimetry method that can be applied to 
microscopic flows and can extract three-dimensional velocities from two-dimensional 
images by using defocused images quantitatively. 

4. Direct visualization of wall-slip in a surface solution 

A microscope-mounted torsional shearing-flow cell was constructed and 
microscopic particle imaging velocimetry employed to directly visualize and map the 
velocity slip layer in a shearing flow of dilute micellar surfactant solutions of 
cetyltriammonium bromide/sodium salicylate. It was shown that the thickness of the 
wall-slip layer is about 1 00 mm at low shear rates, decreasing to around 50-60mm at 
high shear rates. Surprisingly, we found that the wall slip layer emerged only near the 
upper rotating plate of the flow cell, and not the lower surface, a phenomenon that is still 
unexplained. The rheological properties of the wormy micellar solutions were also 
measured by using a stress-controlled rheometer and the results compared to the 
visualization results of wall-slip velocity in the shearing-flow cell. Above the critical 
point, the shear viscosity reaches a plateau and decreases at a higher shear rate where we 



observed that the flow becomes unstable. As the gap of parallel plate geometry or the 
cone angle of cone-plate geometry increased, the critical shear rate decreased, because 
the wall-slip layer formed in the flow cell even below the shear thickening transition. We 
extracted the wall-slip velocities from the rheological data using a Mooney analysis and 
show, apparently for the first time, that the slip velocity from the Mooney analysis is 
consistent with that obtained by direct visualization experiments. 

Key Result : First direct verification that slip velocities extracted from macroscopic gap- 
dependent rheology measurements are consistent with direct microscopic measurements 
of slip, thus validating a 50-year-old technique, the Mooney analysis. Also, microscopic 
imaging of thick slip layers in flowing surfactant solutions. 

Publications 

Four papers have been generated through these studies. Two of them have been 
published and the other two are in preparation for submission. Attachments are the 
published papers and a preliminary, partial, report on micro-fluid dynamics. The latter 
does not yet contain a full write-up of the recently completed FEM study. The follow is a 
list of the papers arising from this grant: 

• Hua Hu and R. G. Larson, Evaporation of a Sessile Droplet on a Substrate, J. 
Phys. Chem. B, 106, 1334-1344(2002). 

• Hua Hu and R. G. Larson, J.J Magda, Measurement of wall-slip-layer rheology in 
shear-thickening wormy micelle solutions, J. Rheol. 46—4, 1001-1021(2002). 

• Hu, H., Larson. R. G., Micro-Fluid Dynamics in an Evaporating Sessile Droplet, 
to be submitted. 

• Chopra, M., Li, L., Hu, H., Bums, M.A., Larson, R.G. DNA Molecular 
Configurations in an Evaporating Droplet Near a Glass Surface, to be submitted 
( 2002 ). 
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The evaporation of a sessile droplet with a pinned contact line is investigated experimentally, by analytic 
theory and by computation using the finite element method (FEM). Because of the low value of R 2 IDt f = 
Cv (l - H)fp = 14 x 10 _5 , where R is the contact-line radius, D is the water vapor diffusivity, c v is the 
saturated water vapor concentration, H is the relative humidity, and p is the liquid water density, the evaporation 
can be considered as a quasi -steady-state process. Hence, the vapor concentration distribution above the droplet 
satisfies the Laplace equation but with a time-varying droplet surface. It is found both theoretically and 
experimentally that the net evaporation rate from the droplet remains almost constant with time for a small 
initial contact angle (0 < 40°), even though the evaporation flux becomes more strongly singular at the edge 
of the droplet as the contact angle decreases during evaporation. We also measured the critical contact angle 
at which the contact line starts to recede and found that it is about 2-4° for clean water on glass. Finally, we 
compare the results obtained by our FEM analysis with an analytical solution and derive a very simple 
approximate evaporation rate expression m(f) = —JiRD( 1 - //)c v (O.270 2 + 1.30), which agrees with the 
theoretical results presented by Lebedev [Lebedev, N. N. Special Functions and Their Application ; Prentice 
Hall: Englewood Cliffs, New Jersey, 1965 and Picknett and Bexon [Picknett, R. G.; Bexon, R. J. Colloid 
Interface Sci. 1977, 61, 366] for any initial contact angle 0 between 0 and jt/ 2 with 6 in radians. The 
approximate expression is also compared with droplet evaporation data from the literature, and good agreement 
is found without any parameter fitting. 


1. Introduction 

The evaporation of a sessile droplet is not only important in 
many heat transfer applications but is associated with common, 
everyday phenomena, such as the annoying ring-like spots left 
on dishes that are allowed to dry. Recently, important new 
applications of this simple phenomenon have emerged. Jing and 
co-workers 1 have developed a high-throughput automatic DNA 
mapping method based on drying droplet. In this technique, 
water evaporation is used both to induce a microscopic flow 
that stretches DNA molecules and to deposit those molecules 
onto a substrate where they can be subjected to a restriction 
digestion. The locations of digestion sites along the DNA strand 
observed in an optical microscope then constitutes on “optical 
map” of the DNA molecule. Droplet drying is also important 
in the creation of arrays of DNA spots for gene expression 
analysis. Li et al. 2 found that the DNA’s stretching behavior is 
strongly affected by the evaporation rate of droplet. At a low 
evaporation rate, the DNA molecules are less stretched and their 
molecular conformations are folded or coiled, whereas at a high 
evaporation rate, DNA molecules are more stretched and their 
shapes become dumbbells or half-dumbbells. Thus, there are 
new motivations for revisiting the old problem of a drying sessile 
droplet. 

The evaporation of a sessile droplet has been studied by Birdi 
and Winter, 3 who obtained the evaporation rate by measuring 
the change of weight of droplets of water on a glass surface 
and concluded that for most of the lime during the evaporation 
the evaporation rate remains constant. Later, they 4 reported that 
the rate of evaporation of sessile drops of w'ater on glass (contact 
angle 6 = 41°) and n-octane on Teflon surfaces is constant in 
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time with the contact line pinned. Shanahan and Bourges 5 - 6 
investigated the evaporation of droplets of water from smooth 
polyethylene and from both smooth and rough epoxy resin 
surfaces. They measured the change of drop height, contact 
angle, and contact-line radius with time and observed that 
evaporation occurs in several distinct stages. In the longest of 
these, the droplet flattens as it evaporates with the contact line 
pinned. Rowan et al. 7,8 presented detailed measurements of the 
change in contact angle and height with time for water on poly- 
(methyl methacrylate), PMMA, and for three alcohols resting 
on Teflon. Starting with a large initial contact angle {0 % 80°), 
both the height and the contact angle decreased linearly with 
time. 

Lebedev, 9 and later Picknett and Bexon, 10 independently, used 
the analogy between diffusive concentration fields and electro- 
static potential fields (they both satisfy Laplace’s equation) to 
the problem of evaporation of a sessile droplet. For diffusion- 
controlled evaporation, the vapor concentration field is equiva- 
lent to the electrostatic potential field around the top half of an 
equiconvex lens. These theoretical results should be valid as 
long as the sessile droplet remains in the shape of a spherical 
cap. Bourges and Shanahan 6 proposed an evaporation model 
for the sessile droplet by taking the concentration gradient to 
be that for a hemispherical droplet of same radius as that of the 
sessile droplet. This approximation is not accurate for a flat 
droplet because the distribution of the evaporation flux along a 
sessile droplet surface is not uniform as it would be for a 
hemispherical droplet. Rowan et al. 7 analyzed the problem 
theoretically using a vapor-phase diffusion model suggested by 
Birdi et al. 3 and derived an approximate analytic equation for 
the evaporation rate. Their model Fits experimental results very 
well for droplets with large initial contact angle. Erbil et al. 11 
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modeled the droplet as an ellipsoidal cap, defined by three 
parameters. By adjusting these parameters, they obtained good 
agreement with experiments. Meric and Erbil 12 reported another 
model considering a “pseudo-sphcrical-cap” geometry (for 
which the droplet height is given by /i 0 = olR tan 61 2, where R 
is the droplet radius and a is an adjustable flatness parameter, 
with a = 1 for a spherical cap), which they argued provided 
much better fits to experimental results. Deegan et al. 13 ' 14 
presented an analytical solution for a droplet with the shape of 
a spherical cap. Although authors neglected the evaporation flux 
distribution along the droplet surface, Deegan et al. 14 used the 
exact analytic expression for the evaporation flux distribution 
derived by Lebedev. 9 However, Deegan et al. did not study the 
relationship between the evaporation rate or flux distribution 
and the contact angle. The purpose of this article is to develop 
a simple, yet accurate, model for evaporation of a small droplet 
with the shape of a spherical cap and to lay the basis for 
calculating the Marangoni force induced by a nonuniform 
distribution of the evaporation flux and for developing an 
accurate model for the complete flow field in an evaporating 
droplet, to be presented in a forthcoming paper. 15 This flow 
model will, in future work, be used for predicting stretching 
and deposition of DNA molecules in an evaporating droplet. 

In this paper, we first use a finite element method to solve 
for the outer vapor concentration and evaporation flux. The 
evaporation flux distribution along the droplet surface is not 
uniform when a droplet is placed on the surface. The nonuniform 
evaporation flux eventually affects the droplet surface temper- 
ature distribution and therefore generates a surface tension 
gradient along the droplet surface, which may be a possible 
contributing factor for the contact-line pinning during evapora- 
tion. Second, we describe a particle-tracking method to measure 
the time-dependent change in droplet shape and the rate of 
droplet evaporation. Then the experimental results are presented 
and compared with the results computed by the FEM method. 
Finally, our model for the rate of the evaporation is compared 
to the results reported in the literature, and a simple, yet accurate, 
empirical expression for the evaporation rate as a function of 
contact angle is obtained. The empirical expression is also 
compared to the theoretical results derived by Lebedev 9 and 
Picknett and Bexon. 10 


z 



Figure 1. Droplet with the shape of a spherical cap rests on a flat 
surface. The contact angle is 0, the local height is h{r,t) y and the local 
evaporation flux is 

of 0.8- 1.0 mm and heights of about 0.3 mm and for slow flows 
(around 1 pm/sec), the Bond number is in a range of 0.03- 
0.04 and the capillary number is around 10~ 8 , so that the droplet 
shape can be regarded as a spherical cap. 

In Figure 1, a small drop of water on a glass surface, whose 
shape is that of a spherical cap, is presented. A cylindrical 
coordinate system is used with radial coordinate r and axial 
coordinate z. S = {h(r,t)\r < R} defines the surface of the 
droplet, where h(r,t) is 

h(r y t) = V /f 2 /sin 2 * 6 — r 2 — tf/tan(0) ( 1 ) 

where 9 is the contact angle and R is the contact-line radius. 
The volume of the droplet is 


jtJKO.OPU 2 + /i 2 (0.r)] 
V(0=— 


( 2 ) 


where /i(0,f) is the droplet height as a function of time and R is 
the contact-line radius. 

From eq 1, when we let r — 0, we have /i(0,f) 

h(0,t) = R tan[0(/)/2] (3) 

Because water evaporates into the ambient air, the vapor 
concentration is distributed nonuniformly above the droplet. At 
the interface between the liquid and the vapor, the vapor 
concentration c is assumed to equal the saturation value c y . Far 
above the droplet, the vapor concentration approaches an 
ambient value Hc y (where H is the relative humidity of the 
ambient air). The difference in water vapor concentration 
c v (l — H) drives the evaporation of water into the air, according 
to the diffusion equation 


2. Theory 

From our experimental observations, we have found that, 
when the contact angle is less than 90°, droplet evaporation 
generally includes two main phases. In the first phase, the 
contact angle decreases while the contact line is pinned. In the 

second phase, the contact line recedes while the contact angle 
remains very small. Because the first phase occupies the 90- 
95% of the total drying time, we only consider this phase in 
which the contact line is pinned. In this section, we develop a 
mathematical model and a corresponding FEM solution for the 
droplet evaporation rate. 

2.1. Mathematical Model. Here, we consider a sessile droplet 
having the shape of a spherical cap resting on a flat substrate. 
The droplet shape is controlled by the Bond number, Bo — 
pgRhJo , which accounts for the balance of surface tension and 
gravitational force on the droplet shape, and the capillary number 
Ca = pu/a, which is the ratio of viscous to capillary forces. 
Here, p is the fluid density, g is the gravitational constant, R is 
the contact-line radius, h 0 is the initial height of the droplet, o 
is the air— water surface tension, p is the liquid viscosity, and 
u r is the average radial velocity induced by droplet evaporation. 
In our experiments, with small droplets with contact-line radii 


^ = DAc (4) 

ot 

where c is the local water vapor mass concentration and D is 
the vapor diffusivity. The boundary conditions are 

1. r < R, z = h(r): c = c v 
2. r > R, z = 0: J = 0 
3. r = z = °°: c~ // v c v (5) 

Here J is the vapor mass flux that is due to evaporation. 

The time required for the vapor-phase water concentration 
to adjust to changes in the droplet shape is of the order of R 2 / 
D, where D is the diffusivity of the vapor in air and R is the 
contact-line radius. The ratio of this time to the droplet 
evaporation time tf is R 2 fDu % c v (l — tT)!p. In our experiments, 
we can take // = 0.4, c\ = 2.32 x 10~ 5 g/cm 3 , and p = 1 g/cm\ 
so that we obtain R z /Dtf ^ c v (l - H)fp = 0.000 014 « 1. 
Hence, the water vapor concentration adjusts rapidly compared 
to the time required for droplet evaporation, and the water 
evaporation can be considered to be at a quasi-steady state. We 
therefore neglect the transient term in eq 4 and obtain the 
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Laplace equation for the vapor concentration distribution: 

Ac = 0 (6) 

As the droplet evaporates, the surface of the droplet descends 
toward the substrate. For this moving-boundary problem, we 
present a finite element method to solve eq 6 to calculate the 
vapor distribution above the droplet and the evaporation flux 
along the droplet surface. _ 

At the air-liquid interface, the local evaporation flux J(rj) 
is expressed as 

7(r,/) — DVc (7) 

The evaporation rate 7, over the whole surface is expressed by 
the following equation: 

J t = f r 0*n) dr g (8) 

The above derivation assumes that the evaporation is not so 
rapid as to alter the droplet temperature enough to change the 
values of the parameters c v or D. 

2.2. Finite Element Method. We apply the FEM method to 
calculate the vapor concentration distribution above the droplet. 
For simplification, we describe all of the detailed derivations 
of FEM model in the Appendix. The final formulations of the 
FEM model are expressed by eqs A-5 to A- 12. Incorporating 
the boundary conditions listed in eq 5, eqs A-5 to A- 12 can be 
numerically solved to obtain the vapor concentration distribution. 

3. Experimental Method 

We use a microscopic particle tracer method to measure the 
rate of droplet evaporation. Figure 2 is a schematic of the 
imaging system used in our experiments to measure the droplet 
evaporation, with Figure 2b as an enlargement of element A in 
Figure 2a. A Nikon inverted fluorescence microscope (Eclipse, 
TE200) with a motorized X-Y stage (Prior, Inc.) is employed 
to locate the fluorescent particles. As shown in Figure 2b, a 
droplet of volume about 0.5 pL is deposited on a clean glass 
cover slip (Dow-Corning, No.l ; 22 mm) and initially sealed by 
a cylindrical cap to block its evaporation (Once the droplet size 
and location have been precisely measured, the cap is removed 
and replaced by a cylinder open at the top to allow evaporation 
to begin while suppressing the effect of air currents on 
evaporation.) The glass cover slip with its holder is placed on 
the motorized X-Y stage. A 40 x objective is coated with water, 
and the water coating is brought into contact with the cover 
slip, thus eliminating the air gap between the sample and the 
objective and to minimize the effect of refraction. Fluorescent 
particles 0.75 /tm in diameter (Polyscience, Inc.) are used as 
tracers to map the droplet profiles at different times and thereby 
calculate the residual droplet volume. Particle images are 
detected by a CCD camera (PC-26C, Super Circuit Co.) with a 
resolution of 640 x 480 pixel Then, the images are recorded 
onto a computer disk and analyzed by the SimplePCI image 
processing system (Compix, Inc.). 

With the cylindrical cap in place to prevent evaporation, we 
move the X— Y stage around the perimeter of the droplet and 
measure the coordinates of the particles at the edge of the 
droplet. The coordinates of these particles are Fitted by a circle 
so that the center and the radius of the droplet are determined. 
From the experimental results for about 500 droplets, we can 
conclude that the contact line of the droplet is indeed a circle 
and the radius of the contact line is about 850 ± 10 microns. 
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l.Glass coverslip holder 2. Fluorescence objective(40X) 

3. Water 4. Glass coverslip 5. Droplet 6. Sealing cylinder cap 


b 

Figure 2. Experimental setup. Figure 2a is a Nikon fluorescence 
microscope, and Figure 2b is a blow-up of part A in Figure 2a. 

We also obtained the residual droplet volume by measuring 
the droplet surface profiles at different times by finding the 
coordinates of fluorescent particles on the surface of the droplet. 
The results are shown in Figure 3, in which the symbols are 
the positions of fluorescence particles on the droplet surface 
and the line is a fit by a circular arc, representing a spherical 
cap. The standard error of the Fit is about 1-2 micron. This 
implies the droplet shape is not affected by adding the tracer 
particles. The tracer particles might enhance pinning of the 
contact line, but it is hard to examine it. Because the droplet 
remains in the shape of a spherical cap at all times during droplet 
evaporation, we need only measure the height of the droplet as 
a function of time in order to reach an accurate measurement 
of the droplet surface and volume. 

4. Results 

4.1. Convergence of the FEM Computation. The conver- 
gence of the FEM computation is tested by systematically 
refining the mesh at the edge of the droplet according to the 
method described in Appendix A. 4. The mesh near the edge of 
the droplet is refined several times over so that the number of 
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r (pm) 


Figure 3. Measured droplet profiles at different times, for water 
evaporating from a droplet of initial radius R — 0.85 mm and height 
/io = 0.329 mm. The symbols show the locations of fluorescent particles 
on the droplet surface, and the lines are the fittings of circular arcs to 
these data. The fitting errors are about 1—2 /tm. 



Figure 4. Relative error vs the numbei of FEM elements. 

elements becomes large and the size of the elements near the 
edge becomes smaller. Subsequently, the evaporation flux and 
the total evaporation rate are calculated using eqs 7 and 8 for 
the different refined meshes. Numerical calculations show that 
as the number of the elements increases the total evaporation 
rate decreases for small numbers of elements and levels off when 
the number of elements exceeds about 5000. The relative error 
defined in eq A-13 is plotted in Figure 4, which shows that the 
relative error decreases linearly with the number of elements, 
as expected for linear shape functions. Figure 4 show’s that when 
the number of elements in the circular sector is larger than 
10 000 the relative error is small enough that the convergence 
criterion is satisfied. We can conclude that after the mesh has 
been locally refined 6-fold according to the method given in 
Appendix A.4, an acceptable FEM computational accuracy is 
attained. Using this protocol, the vapor concentration distribu- 
tion, the evaporation rate, the height of the droplet, the contact 
angle, and the droplet volumes at different times are computed, 
as described in the following. 

4.2. Vapor Concentration Distribution. When using the 
finite element model, eqs A-5 to A- 12, to solve the vapor 
concentration distribution, we need to consider the boundary 
condition in eq 5, which is c = Hc\ at r,z °°. We impose this 
condition along a boundary at (r 2 + z l ) m = KR , where K is a 
constant much larger than unity. By testing different values of 
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Figure 5. Contour plot of the vapor concentration distribution above 
a droplet of radius R = 1 mm and height ho = 0.364 mm. The 
parameters used in the FEM method are vapor diffusivity D = 26.1 
mm 2 /s, relative humidity H — 0.40 (i.e., 40%), and saturated vapor 
concentration on the droplet surface r v = 2.32 x 10~ 8 g/mm\ which 
is the value obtained form the CRC Handbook 16 at 25 °C. These 
parameter values are also used in Figures 6— 10. The gray bars represent 
the vapor concentration in g/mml 



Figure 6. Vapor concentration distribution along the directions 2 = 0 
and r = 0 from FEM (solid lines) and analytic solution (dashed lines). 
The FEM results along each direction superimpose on the analytic 
results. 


K between 10 and 100, we find that the calculated vapor 
concentration at position (r 2 + z 2 ) l/2 = 20 R is close enough to 
the ambient vapor concentration c«, that the deviations are 
negligible; that is, (c — c«,)/(c v “ c\») < € % 0.002. Therefore, 
in our FEM analysis, we use the boundary (r 2 + z 2 ) m = 20 R 
to approximately represent r,z — * °°. 

The vapor concentration distributions at different times are 
then computed by the FEM method and presented in Figure 5. 
Figure 5 is a contour plot of the vapor concentration distribution 
when the droplet is starting to evaporate. The contour lines are 
concentrated near the droplet surface, where the vapor concen- 
tration has a large gradient. The vapor concentration along the 
lines r — 0 and z = 0 is plotted in Figure 6. When r and z are 
greater than 5/?, the changes in vapor concentration are small, 





1338 J. Phys. Chem. B, Voi 106, No. 6, 2002 


Hu and Larson 


the evaporation flux along the droplet surface 



Figure 7. Evaporation flux along the droplet surface. The insert shows the magnitude of the evaporation flux along the droplet surface. The flux 
vector at the edge of the droplet does not orient along the normal direction, because when the mesh is refined several times near the edge of the 
droplet, the edge of the last element at r — R, which has become very small, does not perfectly match the droplet surface profile. This artifact is, 
however, confined to a single element and has negligible effect on the overall flux. 


and the normalized vapor concentration, (c* - Hc v )/( 1 — //)c v , 
is so small, 0.02, that the cutoff, ( r 2 + z*) m = 20 R, is acceptable. 

4.3. Vapor Flux above the Droplet. From the vapor 
concentration distribution, we calculate the vapor flux from eq 
7. Along the droplet surface, the flux increases as one moves 
from the center top of the droplet to the contact line at the edge, 
where the flux is theoretically infinite, see Figure 7. The total 
evaporation rate can be evaluated from eq 8 by integrating along 
the droplet surface. The evaporation rate for a series of droplet 
shapes is used to compute the change in droplet volume vs time 
in the next section. 

4.4. Droplet Volume, Height, and Contact Angle vs Time. 

By repeating the FEM analysis for a series of droplet heights, 
we simulate the droplet evaporation process. The radius and 
the initial height of the droplet are 1 and 0.364 mm. The 
parameters used in the FEM analysis are the vapor diffusivity 
D = 26.1 mm 2 /s (from CRC Handbook of Chemistry and 
Physics 16 ), the relative humidity H = 0.4, and the saturated 
vapor concentration c v = 2.32 x 10 8 g/mm 3 . A small time 
step of about 0.02/f is used lo calculate the time-dependent 
volume. At each time step, the loss of water is determined from 
the product of the total evaporation rate integrated over the 
droplet surface and the time step. As we described in section 2, 
the evaporation rate can be obtained from the vapor concentra- 
tion distribution by using eqs 7 and 8 on the droplet surface. 
The new droplet volume is thereby calculated from the loss of 
the solvent and the previous droplet volume, and the new droplet 
surface profile can be derived from eq 2, for a spherical cap. 
This procedure is accurate, because the capillary number is low 
(Ca 10~ 8 ; so that the droplet remains a spherical cap), the 
contact line is pinned, and the vapor concentration field is quasi- 
steady; that is, it adjusts rapidly whenever the droplet shape 
changes. Thus, we can convert this moving- free- surface problem 
into a simple series of solutions to Laplace’s equation. The 
height of the droplet and the contact angle decrease roughly 
linearly with time as shown in Figure 8, suggesting a nearly 
constant evaporation rate. The total evaporation rate, shown in 
the insert to Figure 8, is not quite constant; it decreases slightly 
during the droplet evaporation for an initial contact angle of 
40°. 

4.5. Experimental Results. In our experiments, we obtained 
the droplet evaporation rate by measuring droplet height at a 
series of times. Doing so is reliable because, as shown in Figure 
3, the droplet shape remains a spherical cap during evaporation. 
From the height and the contact-line radius, the droplet volume 
at different times is calculated. The residual droplet volume vs 



Figure 8. Height of the droplet and the contact angle vs time. The 
insert shows the evaporation rate vs time. 

time is plotted in Figure 9 and compared with the results of the 
finite element simulation using the experimentally derived 
parameters. We can see that there is very good agreement 
between the experiment and the FEM calculations. The FEM 
results show that the evaporation rate changes slightly during 
evaporation as seen in Figure 8 and the solid line in Figure 9 is 
not a straight line. It changes its slope slightly during the initial 
drying process but has become almost constant near the end of 
drying. 

In the experiments, we also measured the critical contact 
angle, the angle at which the contact line starts to recede for 
about 50 droplets, and obtained an average of 2-4°. 

5. Discussion 

5.1. Approximate Expression for the Evaporation Rate. 
The excellent agreement between the FEM computation and 
the experimental measurements confirms that the droplet 
evaporation is a quasi-steady-state process. This should be true 
whenever R 2 IDt{ ^ c v (l — H)/p & cjp is small; that is, 
whenever the vapor phase has a density much smaller than that 
of the liquid, which is always true except near supercritical 
conditions. Thus, the quasisteady approximation for the vapor 
concentration field should be valid even for very rapidly 
evaporating droplets. It should be kept in mind, however, that 
for very rapidly evaporating droplets large temperature non- 
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Figure 9. Symbols giving the residual droplet volume vs time 
calculated from the height of the droplet at different times and the 
contact-line radius. The thin solid line is the result of the FEM 
calculations, using the experimental conditions as parameters, i.e., R 
= 0.95 mm, h 0 = 0.364 mm, vapor dtffusivity D = 26. 1 mm 2 /s, relative 
humidity H = 0.38 (i.e., 38%), and saturated vapor concentration on 
the droplet surface c\ - 2.32 x 10 8 g/mm 3 at temperature 25 °C. The 
dashed line is calculated by the approximate evaporation rate expression 
eq 21 using the same values of parameters as the FEM method. The 
thin dashed line is calculated by Picknett and Bexon’s model 9 using 
the same value of parameters as the FEM method. 
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Figure 11. Evaporation rate exponent A as a function of the contact 
angle, obtained by fitting eq 9 to the FEM results. 

but starts to decrease linearly with 6 for 0 less than 60°. At 6 
= 0°, the ratio reaches the value of 0.6377. The error between 
the FEM result and the exact solution shown by the dashed 
line in Figure 10 is less than 1.5%. Surprisingly, Figure 1 1 shows 
that A decreases linearly with 0 over the range of 0—90°, and 
so A can be expressed by a simple linear function, X(6) = 0.5 
— 6/tt , which is different from a result cited by Deegan et al., 
X(6) = ( 7 i - 26)f(2jt - 20). However, this latter formula for 
X(6) applies to the equation J(r) ~ (I — r)~ x and not to the 
equation used here (and by Deegan et al.), namely, J(r) (1 - 
r 2 ) - ^. Obviously, the formula for the exponent X(6) will differ 
for these two expressions. By using the formula X(0) = 0.5 — 
Otei , the largest relative error between the prediction from eq 9 
and the results from the FEM analysis over the range of 0-90° 
is less than 6%, and this value decreases toward zero as the 
contact angle approaches either 0° or 90°. From Figures 10 and 
1 1, we conclude that the evaporation rate is a function of the 
contact angle 6. Using eq 9, the evaporation rate becomes 

~m(0 = f s 0(r,tYn) dS - 

fg2jl r-R 2 'l 0 (8)(\ — + 1 dr (10) 



Contact angle (radians) 

Figure 10. Ratio Jd(9)Uo{Jif2) versus the contact angle. The solid line 
is calculated by FEM analysis, whereas the dashed line is from the 
analytical solution. 

uniformities may develop because of latent heat and this will 
affect the evaporation flux. 

The evaporation flux on the droplet surface calculated in our 
FEM analysis can be fitted by the following equation suggested 
by Deegan et al.: ,3 U 

(7-«) = y 0 (i - f 2 )' 1 (9) 


where S is the area of integration and Sh(rj)/dr is the derivative 
of h(r y t) with respect to r, given for a spherical cap by 



We find that the term [(d/i(r, t)/dr) 2 + \\ [a in the integration 
kernel in eq 10 can be approximated by 

JjffiTl = ( 1 - sin 2 6 r 2 )-° 5 * ( 1 - r 2 )"™ ( 1 2) 


where A is a fitting parameter representing the nonuniformity 
of the evaporation flux on the droplet surface and r s r/R. From 
fits of eq 9 to our FEM results, we obtain values of Jo and A 
for different contact angles, see Figures 10 and 11. Because 
Deegan et al. did not precisely define the relationships between 
J 0 , A, and the contact angle 0 , we determine Jq{6 ) and X(6) 
versus the contact angle empirically by using the finite element 
method. From Figure 10, we find that the ratio J 0 (6)/Jq(jt12) is 
close to unity when the contact angle 0 is in a range of 70-90° 


where 6(6) is an empirical function of contact angle 0. Inserting 
eq 12 into eq 10, we have 

—m(t) = rR 2 -J 0 (6)( 1 - r 2 ) _A(#) dr (13) 

where A(6) = X(6) + 6(6) is a combination of two factors, 
one is a parameter reflecting the nonuniformity of the evapora- 
tion flux and the other reflects the droplet surface area per unit 
area of substrate at each value of r. When the contact angle is 
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Contact angle (radians) 

Figure 12. A(0) versus the contact angle 0. The solid symbols are 
the slopes of the curves in the insert obtained by fitting ln[tf(r,6i)] = 
ln{( 1 - r 2 ) _ ^ o> [(a/i(r ? 0)ar) 2 + l] l/2 } with -A(0)ln(l - r 2 ). The line is 
a parabola, eq 14, fitted to the solid points. 

90°, because \{nl 2) is 0 and 6{nf2) is 0.5, A(jt/ 2) equals 0.5. 
When the contact angle is 0°, because A(0) is 0.5 and (5(0) is 0, 
A(0) also equals 0.5. Using the results of our FEM analysis, 
we fit the term (1 — r 2 ) _A<0) [(3/j(r, t)/Sr) 2 + \] ia in eq 10 with 
the equation (1 — r 2 )* A(tf) at a given contact angle to obtain 
A(0) and plot this in Figure 12 as solid diamonds. In the insert 
to Figure 12, we plot the term ln[ = ln{( 1 — r 2 ) 
[(Qh(r,0)/dr) 2 + \] ia } against — ln( 1 - r 2 ), and the slopes of 
the lines give A(0). The points in Figure 12 can be fitted well 
by a parabola, namely, 

A(0) = 0.2239(0 - .t/4) 2 + 0.3619 (14) 

where 6 is given in radians. 

Substituting eq 14 into eq 13, and integrating, gives 

_ 7tR 2 J (i ( 0) _ JI R 2 Jq(0) 

~~ 1 - A (0) 1 - (0.2219(6 - nl\f + 0.3619) 

(15) 

Equation 15 is a general formula tor the droplet evaporation 
rate for any contact angle (0° < 9 < 90°). In eq 15, we note 
that J 0 (0 ), the evaporation flux at the center of the droplet, is 
also a function of the contact angle 0, as shown in Figure 10. 
To simplify eq 15, we calculated the ratio of Jq(9)/Jo(tz/2)(\ — 
A (6)) using the FEM results and plot this in Figure 13, which 
is given by solid squares. In Figure 1 3, the line is obtained from 
the fitting function 

JJ0) 7 

— — — — = JJx/2)(0.278~ + 1.30) (16) 

1 - A (6) 0 

where 0 is given in radians. The results predicted by eq 16 are 
compared to the exact solution for contact angle ranging from 
0 to 90°, and the largest error is less than 6%. 

Substituting eq 16 into eq 15, we derive an approximate 
expression for the droplet evaporation rate at any contact angle 
(0° < 6 < 90°), which is 

—m(t) = nR 2 J Q (nl2 XO.270 2 + 1.30) (17) 

where Jq{7i/2) is the evaporation flux for the contact angle 90°. 


Hu and Larson 



Contact angle 0 (radians) 


Figure 13. Ratio - A (9)) versus the contact angle, 

calculated from Figures 10 and 12. 

F or 0 = 90°, the evaporation flux is uniform everywhere along 
the droplet surface. The solution of eq 6 for a contact angle of 
90° is, by symmetry, the same as the solution for a droplet 
suspended in the air without the presence of a substrate. For a 
suspended droplet, the boundary conditions are (r 2 4* z 2 ) ,/2 = 
R:c = c\ and (r 2 + z 2 ) 1 ' 2 = = Hc v , and we derive the 

evaporation flux Jq{tiI2) from eq 7, which is 

D( 1 - //)c v 

J 0 (7i!2) = (18) 

Combining eqs 17 and 18 gives 

-m(t) = tzRD(\ - W)c v (0.276 2 + 1.30) (19) 

From eq 19, we can see that at a given contact angle the 
evaporation rate is proportional to the contact-line radius R, the 
vapor concentration difference (1 — H)c v > and the diffusivity 
D and depends weakly on the contact angle 6. Figure 13 shows 
that when the contact angle is less than 40° (0.7 rad) the 
dependence on 6 almost becomes flat, and therefore, the 
evaporation rate is almost constant. Equation 19 agrees well 
with the theoretical results obtained by Picknett and Bexon. 10 
We calculate the evaporation rate as a function of the contact 
angle according to eq 19 and compare the result to that derived 
by Picknett and Bexon. The averaged relative error between 
the two predictions is less than 1%. However, eq 19 is different 
from the prediction of the model presented by Bourges and 
Shanahan. 6 When the contact angle 6 is close to 0°, their model 
predicts an evaporation rate of 0, which is incorrect. This is 
because that Bourges and Shanahan assume that the vapor 
concentration gradient is uniform along the surface of the sessile 
droplet and is set by the radius of curvature of the droplet 
surface, which becomes infinite when the droplet becomes flat, 
leading to zero evaporation flux. Actually, from our FEM results, 
Figure 7, we can see the concentration gradient along the surface 
of a sessile droplet is not uniform compared to the uniform 
concentration gradient on the surface of a spherical droplet. 

Equation 19 gives the exact solution for two limiting cases, 
for a contact angle of 0° and a contact angle of 90°. When the 
contact angle is close to 90°, eq 19 gives 

—m(t) = 2 jtD(\ — H)c xr R 


( 20 ) 
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whereas when the contact angle 0 ss close to 0°, eq 19 gives 

-m(t) = AD(\ — H)c y R (21) 

Both eqs 20 and 2 1 are consistent with the results of Picknett 
and Bexon 10 for the two cases of the contact angles, 90 and 0°, 
respectively. 

We use eq 21 to calculate the dioplet volume vs time and 
compare this with the results from the FEM method and the 
experimental measurement as shown in Figure 9. Figure 9 shows 
that the results calculated by eq 21 are consistent with those 
from the FEM method, Picknett and Bexon’s model, and the 
experiments. The relative error between eq 21 and the experi- 
ments is 3.6%. This implies that when the initial contact angle 
is less than 40° the evaporation rate can be approximated by eq 
21 and the evaporation rate can be regarded as a constant. The 
prediction from the FEM analysis, which is the thin solid line, 
is very close to the prediction of Picknett and Bexon’s model, 
which is the thin dashed line. The relative error in the drying 
times predicted by FEM and the theoretical model of Picknett 
and Bexon is less than 1.3%. 

5.2. Comparison with the Exact Solution. Deegan et al. u 
reported an analytical solution for the Laplace equation, which 
was derived by Lebedev. 9 Lebedev considered a charged surface 
formed by the union of two spherical domains, and derived the 
electrostatic potential distribution in the space exterior to this 
surface by using toroidal coordinates. In toroidal coordinates, 
the potential distribution is given by 


droplet. We consider the electrostatic potential u as a dimen- 
sionless vapor concentration, u = (c ~ t\x,)/(t’ v — c«), where c\ 
and are the vapor concentration on the droplet surface and 
in the ambient, respectively. Then u = l on the droplet surface. 
For the lens-like geometry, we let P\ = n ~ 0 and pi = n + 
0 , where 0 is the contact angle. From eq 25, we obtain the vapor 
concentration distribution 



y/2 cosh a — 2 cos pj^ 


cosh(0r) cosh[(2jr — p)r] 
cosh(^r) cosh[(jr — 0) r] 


^-u/W cosh a > dr < 26 ) 


Once the vapor concentration is known, from eq 7, we can 
obtain the evaporation flux distribution along the droplet surface, 
which in toroidal coordinates is 


- _ D(cosh a — cos p) 

(J*n) — — 


(27) 


Inserting eq 26 into eq 27, the evaporation flux is 


(J-n ) : 


D(c v - cJ 

sin 0 

R 

2 


- + (cosh a + 


cos 


0f n) J, 7 ~ e ) r ! T/, -(i/ 2 )+ir(c°sh a) dr 


cosh(jr r) 


(28) 


u = Vy/2 cosh a - 2 cos pj*^ 


cosh[(;r - £,)r] sinh[{^ - fi 2 ) r] + cosh[( :r - fi 2 )i] sinh[(2;r + /? t - ft) r] 
cosh(^r) sinh[<2;i 4- /?, - fi 2 )t] 

P_ ( ,/ 2 )+ rt ( cosh a) dr (22) 


where a and p are the toroidal coordinates, V is the potential 
on the surface, and pi and pi are two angles of n - 6 and n + 
0, P_ (1/2)+iI <cosh a) is the Legendre function of the first kind 
and is expressed by 

2 sin(r t) 

P-am+i r(c°sh a) = -cothfjnO/ -====== d; 

71 V 2 cosh t — 2 cosh a 

(23) 


which corrects eq A2 of Deegan et al. u 

Eqs 6 and 28 are the expressions for the vapor concentration 
distribution above the droplet and the evaporation flux along 
the droplet surface, respectively. Both of them have complicated 
integrations, so closed forms are not available. A numerical 
method is therefore used to solve these two formulas. 

Now we evaluate eqs 26 and 28 for two special cases, 
0 = 0 and 90°. For the case of 0 = 90°, eq 26 becomes 


c 



4 -, 


cosh a — cos P __ R 
cosh a -f cos p r 


(29) 


where r = (r 2 + z 1 ) if2 
Equation 29 becomes 


Although eq 22 has a different form from Picknett and 
Bexon’s solution, they are identical because they are derived 
from the same model. 

The toroidal coordinates a and P are related to the cylindrical 
coordinates r and z by 


R sinh a 

(24) 

cosh a — cos P 


R sin P 

(25) 

cosh a — cos p 



where R is the contact-line radius. 

The electrostatic potential around a lens-shaped object with 
uniform surface potential and the vapor concentration distribu- 
tions above an evaporating droplet are equivalent, because they 
are both solutions of the Laplace equation. By symmetry, the 
solution to Laplace’s equation in the half-plane above a spherical 
cap is identical to that of a lens-like shape formed from two 
back-to-back spherical caps. Thus, we can apply the solution 
of the electrostatic potential in toroidal coordinates for a lens- 
like shape to the vapor concentraiion distribution above the 


(J-n) 



(30) 


Equations 9 and 30 are consistent with the results obtained by 
solving the Laplace equation in cylindrical coordinates. Inserting 
eq 30 into eq 10, the evaporation rate obtained is identical to 
eq 20. 

For the case of 6 = 0°, eqs 26 and 28 can be solved: 


c - 3 i 

= - H — arctan 

c v — c' 2 it 


cos p — sinh 2 |^j 


cos|^jV2 cosh a - 2 cos p 


(31) 


D ( c v - O 2 

(J-n) = 


R 


-cosh(^ 

JT \2j 


In cylindrical coordinates, eq 32 becomes 

- D(c * " cJ 2 n r-V 05 

( J -n) -<l-r) 


(32) 


( 33 ) 
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Equation 33 implies that the fitting parameter X in eq 9 
becomes 0.5 when the contact angle reaches the limiting value 
of 0°. From eqs 30 and 33, we find the ratio Jo(0)//o(^/2) = 
[(2fji)(D(c\ - c„)/R)]1(D(l\ - c~)/R) = 2 /ji. The ratio from 
Figure 10 computed by the FEM method for 0 — 0 is 0.6377, 
which is indeed very close to the value 2 hi = 0.6366. 

Substituting eq 33 into eq 10, we also rederive the same 
evaporation rate for 0 — 0° given by eq 21. 

We can now compare the FEM results with the analytical 
results obtained by solving eqs 26 and 27. In Figure 6, the vapor 
concentration distributions along the r = 0 and z — 0 directions 
calculated by FEM are indistinguishable from the results 
obtained from eq 26. Figure 10 shows that the ratio Jo(0)fJ<y 
(ji/2) computed by FEM and from the analytical solution are 
nearly identical. 

We note that the analytic solution is only available for the 
special case of a spherical cap, which is valid for small droplets 
(R < 1 mm). For larger droplets, the bond number exceeds 0.1, 
and gravitational sag becomes important. The analytic solution 
also fails when the evaporation rate is fast enough to change 
the temperature of the droplet enough to affect the saturated 
vapor concentration cv. In the case of gravitational sag, eq 1 
for the droplet surface profile must be replaced by an ordinary 
differential equation for the static droplet shape under the 
influence of gravity and surface tension, which must be solved 
each time step. If temperature variations become large, the 
temperature field inside the droplet must be obtained each time 
step by solving a quasisteady heat equation (as is done in a 
forthcoming paper 15 ). In both cases, however, the most important 
simplifications of our approach remain valid: the vapor 

concentration field is at quasi -steady-state (because R 2 /Dt{ = 
c v ( 1 - H)!p is small), and the droplet shape is at static 
equilibrium (because Ca = //u Jo is small). Thus, the simple 
quasi-steady-state FEM method presented here remains accurate 
for a very wide range of conditions, including conditions for 
which eq 26, the electrostatic solution, does not apply, such as 
when the droplet is not a spherical cap (because of gravity) or 
the temperature is not uniform within the droplet because of 
rapid evaporation. 

5.3. Comparison of Model Predictions with Other Ex- 
perimental Results. We now compare the results calculated 
from the simple empirical formula eq 19 (which fits the FEM 
and analytic results almost perfectly) with the experimental 
results reported by Birdi et al7 and Rowan et al. 7 Birdi et al. 3 
studied the evaporation of droplets of water on glass by weighing 
the residual droplet mass. In general, we should use eq 19 to 
compare with their experiments, but this equation is well 
approximated by eq 21 when the initial contact angle 6 is les 
than 40°, as is the case in these experiments. Their paper does 
not give the parameters D , //, and t v . Therefore, we obtain the 
experimental term 4D(1 - H)c v from one of their experiments 
and then apply this constant to calculate the results for other 
droplets with different contact line radii. The results are plotted 
in Figure 14, which shows that if we arbitrarily fit the set of 
experimental data for the droplet with contact-line radius of 2.01 
mm by eq 21 we obtain the fitting constant 4D(1 - H)c v = 6.5 
x 10~ 5 g s" 1 mm -1 with correlation coefficient R = 0.9999. 
With this constant, the residual masses of the droplets of radii 
2.53 and 2.93 mm are calculated by the integration of eq 21 
with time. The results are plotted in Figure 14, in which the 
solid lines are the calculated results, which are very close to 
the experimental results. The average relative errors between 
the theoretical predictions and the experiments are 3.1% and 
3.5% for droplet radii 2.53 and 2.93 mm, respectively. 



Figure 14. Comparison of the time-dependent weight from the model 
predictions (eq 21) and the results published by Birdi et al7 for droplets 
of water of radii 2.01, 2.53, and 2.93 mm on glass at T = 22 °C. The 
theory for R = 2.01 mm was fit to the data by adjusting 4D(l - //>\ 
to the best-fit value of 0.000 065 g mm" 1 s' 1 , and this was held fixed 
for the other experiments. The symbols are experimental results and 
the lines are the model predictions. 



Figure 15. Comparison between the model prediction (eq 19) and the 
results published by Rowan et al. 7 for water droplets with radii R ~ 
0.585, 0.491, 0.451, 0.381, 0.324, and 0.293 mm on PMMA substrates 
at T = 21.5 °C. All parameters in eq 19 are obtained from Rowan’s 
paper, and they are vapor-phase water diffusivity D— 17 mm 2 /s, relative 
humidity H = 0.55 (i.e., 55%), and saturated vapor concentration on 
the droplet surface c y = 1.9 x 10" 8 g/mm 3 , The symbols are the 
experimental results, and the lines are the model predictions. 

We also compare our model's predictions with Rowan’s 7 
results in Figure 15. Rowan et al. studied the height and the 
contact angle as functions of time for a water droplet on a 
PMMA substrate by using a Kriiss contact-angle meter. In their 
experiments, the droplets have an initial contact angle of about 
80°. Therefore, we employ eq 19 to calculate the residual droplet 
volume versus time. Here, the parameters in eq 19 are directly 
obtained form Rowan’s paper. They reported that the water 
vapor diffusivity is D = 17 mm 2 /s (which is a fitting value by 
using their model to their experimental results), the relative 
humidity is H = 0.55 (i.e. 55%), and the saturated vapor 
concentration on the droplet surface is c\ = 1.9 x 10 8 g/mm 3 . 
The contact line radii for six different droplets are 0.585, 0.491, 
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0.451, 0.381, 0.324, and 0.293 mm From Figure 15, we can 
see that, when droplet is bigger ( R = 0.585 and 0.491 mm), 
our model predicts the experimental results fairly well, and the 
average relative errors between the predictions and the experi- 
ments are 7.4% and 9.7%, respectively. As the droplet size 
becomes smaller, the relative errors between the model predic- 
tion and the experiments become larger and are in the range of 
10-25%. It is nearly impossible that evaporation cooling 
produces such big errors because in the forthcoming paper 15 
we have calculated that the temperature drops in the droplet 
are only about 0.02 °C, which has hardly any effects on c\ and 
D. Possible reasons for the errors are that when droplet size is 
smaller (below 0.451 mm) the contact angle of the droplet is 
hard to determine precisely by the method described in the paper 
or that the humidity used here is not the true value. 

In all of these cases, the overall evaporation rate remains 
nearly constant as the droplet dries out, and thus, the singularity 
in evaporation flux at the droplet edge is of little consequence 
for the overall rate of droplet drying. However, this singularity 
in evaporation flux creates a droplet llow field with a singularity 
at the contact line. In addition, small temperature variations 
along the droplet surface, which have negligible effects on the 
drying rate, induce Marangoni stresses (Scriven 17 ), which affect 
qualitatively the flow inside the droplet. In our forthcoming 
paper, 15 this flow field will be examined both theoretically and 
experimentally. 

6. Conclusion 

An FEM model is developed to solve the vapor concentration 
distribution and the evaporation flux above a droplet that is small 
enough that its shape is not influenced by gravity and is therefore 
a spherical cap. The vapor phase water concentration field 
adjusts rapidly to changes in droplet height and can be regarded 
as quasisteady. The evaporation flux along the droplet surface 
is not uniform and becomes singular at the edge of the droplet 
and can be fitted by the expression /o(#)(l “ ( rlR) 2 y X(0) 
suggested by Deegan and co-workers, where r is the radial 
position along the droplet, R is the droplet radius, and Jq(9) 
and X(0) are empirical functions of the contact angle 6 , which 
give the precise predictions comparing to the results obtained 
by Deegan and co-workers. 14 Neglecting the nonuniformity of 
the evaporation flux along droplet surface will lead an inaccurate 
theoretical result. FEM results also show that the overall 
evaporation rate is almost constant over the whole evaporation 
period when the initial contact angle is less than 40°. The FEM 
results agree well with the experimental measurements. Finally, 
for contact angles between 0 and 90°, an accurate approximate 
expression for the evaporation rate is presented: — m(/) = 
—jzRD( 1 - H)c v (0.216 2 + 1.30), where D is the gas-phase 
diffusivity, H is the relative humidity, 6 is the contact angle in 
radians, and m(t) is the time-dependent droplet mass. This 
equation is confirmed by Picknett and Bexon’s results. For 0 
< 40°, this expression can be reduced to just — m(t) = 4RD(\ 
- H)C\. The results predicted by this model compare very well 
with literature data without parameter fitting. The FEM results 
are also confirmed by an analytical formula derived by Lebedev 9 
for the electrostatic potential produced by a charged lens and 
by the theoretical solution of Picknett and Bexon. 10 
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Appendix 

We here derive the FEM element model that is used to 
calculate the vapor concentration distribution. 
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A.l. Weak Formulation. From a variational analysis, the 
weak form of eq 6 in the cylindrical coordinate system is 

/ a w(V-(DrVc)) dQ = 0 (A-l) 

where Q is the domain of the vapor phase, which is the 
semiinfinite half plane above the droplet excluding the droplet 
volume, and w is the weighting function. Integrating eq A-l by 
parts gives 

f Q Vw-(DrVc) dQ - f r (rDVc)-n dr = 0 (A-2) 

The second term in eq A-2 is an integration over the boundary 
T, where T = T g + T h , and T g D T h = 0, and n is the outward- 
pointing unit vector along the surface. F g is the droplet free 
surface, and Th is the dry substrate surface excluding the area 
on which the droplet rests. Boundary condition 1 in eq 5 applies 
along T g , and boundary condition 2 applies along Th The 
properties of the weighting function w are as follows: along 
T = r g , w = 0; along T = T/,, w = 1 , which we have used to 
derive eq A-2. 

A.2. Galerkin Method. To complete the finite element 
model, we apply the Galerkin method to approximate the 
concentration c and the weighting function w: 


n 


c= 

A = 1 

(A-3) 

n 

A = 1 

(A-4) 

d A and ** are the coefficients on node A, n is the total number 
of nodes, and Na is shape function for node A. 

When eqs A-2, A-3, and A-4 are combined, the finite element 
model discretization of eq 6 is obtained as follow: 

KD = F 

(A-5) 

In this equation, K, F, and D are 


K=[k AB \ = [f Q WN A -VN B rdQ] 

(A-6) 

F=irj = [/ r W V Hrdr] 

(A-7) 

D = [d A \ = [d v d v d v d A f 

(A-8) 

A.3. Shape Function. We employ linear triangular elements 
in the finite element model, eqs A-6 to A-8. The shape functions 

are 


N t = L t 

(A-9) 

.? 

II 

,r 

(A- 10) 

ii 

af 

(A-l 1) 

L, +£, + £,= 1 

(A-12) 


Lu Li, and L 3 are the local variables in each triangular element, 
and only two of them are independent variables, as specified 
by eq A- 12. 

A.4. Convergence Criteria and Mesh Refinement. Because 
two kinds of boundary conditions meet at the edge of the droplet, 
the evaporation flux becomes singular there. To overcome this 
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difficulty, we refine the mesh near the edge of the droplet. For 
convenience, we use a commercial software package, Ansys 
5.6, to generate the mesh for our FEM analysis. Initially, the 
mesh is generated coarsely, and the FEM code is then used to 
calculate the vapor concentration and the evaporation rate. Next, 
we choose a circular sector with a radius of r\ — 0.5 R centered 
at the edge of the droplet in the initial mesh and refine all of 
the elements in this area. The vapor concentration and the 
evaporation rate are calculated again for the new mesh. We 
compare the current total evaporation rate J u on mesh refinement 
/ to J lti - j, which is for the less refined mesh i— 1 , to see whether 
they satisfy the following criteria: 

e = — < 0.005 (A- 13) 

where J x is expressed by eq 8 

j, = /(?•») dr, (A- 14) 

To obtain an accurate vapor concentration distribution, the 
mesh near the edge of the droplet is refined continuously until 
the criterion in eq A- 13 is satisfied. However, the new, ith, 
refinement is concentrated in a circular sector with radius r, = 
0.68r,-i, where n~\ is the radius of the circular sector used in 


the n-ith refinement. The final refined mesh, after six iterated 
refinements, has 11071 elements in the circular sector. 
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